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Abstract 

Mechanical unfolding of RNA structures, ranging from hairpins to ribozymes, using laser optical 
tweezer (LOT) experiments have begun to reveal the features of the energy landscape that cannot 
be easily explored using conventional experiments. Upon application of constant force (/), RNA 
hairpins undergo cooperative transitions from folded to unfolded states whereas subdomains of 
ribozymes unravel one at a time. Here, we use a self-organized polymer (SOP) model and Brownian 
dynamics simulations to probe mechanical unfolding at constant force and constant-loading rate 
of four RNA structures of varying complexity. For simple hairpins, such as P5GA, application of 
constant force or constant loading rate results in bistable cooperative transitions between folded 
and unfolded states without populating any intermediates. The transition state location (Ax^'^) 
changes dramatically as the loading rate is varied. At loading rates that are comparable to that 
used in LOT experiments the hairpin is plastic with Ax"^^ being midway between folded and 
unfolded states whereas at high loading rates Ax|^'^ moves close to the folded state i.e, RNA 
is brittle. For the 29 nucleotide TAR RNA with the three nucleotide bulge unfolding occurs 
in a nearly two state manner with an occasional pause in a high free energy metastable state. 
Forced unfolding of the 55 nucleotide of the Hepatitis IRES domain Ha, which has a distorted 
L-shaped structure, results in well populated stable intermediates. The most stable force-stabilized 
intermediate represents straightening of the L-shaped structure. For these structures the unfolding 
pathways can be predicted using the contact map of the native structures. Unfolding of a RNA motif 
with internal multiloop, namely, the 109 nucleotide prohead RNA which is part of the 4'29 DNA 
packaging motor, at constant value of occurs with three distinct rips that represent unraveling of 
the paired helices. The rips represent kinetic barriers to unfolding. Our work shows (i) the response 
of RNA to force is largely determined by the native structure; (ii) only by probing mechanical 
unfolding over a wide range of forces can the underlying energy landscape be fully explored. 
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I. INTRODUCTION 



The discovery of self-splicing catalytic activity of ciliate Tetrahymena thermophila ri- 
bozyme 0, 0] and subsequent findings that RNA molecules play an active role as enzymes 

in many cellular processes have revolutionized RNA research. Just as for protein folding 
the structures and functions of RNA enzymes (ribozymes) are linked. As a result, the RNA 
folding problem, namely, how a nucleotide sequence folds to the native state conformation is 
important in molecular biology. Several studies from a number of groups j^ilESSIESl have 
shown that, under in vitro conditions, ribozymes have rugged energy landscapes. Despite 
significant advances in our understanding of how ribozymes fold, several outstanding issues 
remain which can be addressed using single molecule experiments 0, IE IE B B 113 • 

Thermodynamic and kinetic measurements in ensemble and single molecule florescent 
energy transfer experiments are typically made by varying the concentration of counterions. 
Recently, using laser optical tweezer setup mechanical force has been used to trigger folding 
and unfolding of RNA molecules at a single molecule level 3|. Mechanical force, applied 
to a specific position of the molecule, induces sequence and structure-dependent response, 
which is reflected in the force-extension curve (FEC) that is usually fit using the worm-like 
(WLC) model [11111. The stabihty of RNAs is inferred by integrating the FECs. For 



simple motifs, such as hairpins, it has been shown that the stability of the native structures 
can be accurately measured using mechanical unfolding trajectories which exhibit multiple 
transitions between the folded and the unfolded state when the force is held constant (ill ]. 
Similarly, thermodynamics of ribozymes can also be obtained using the non-equilibrium work 
theorem [lE, Il3 • 

Mechanical force has also been used to probe unfolding and refolding kinetics of RNA. 
The cooperative reversible folding of hairpins has been shown by monitoring the end-to-end 
distance (R), a variable conjugate to the mechanical /, as a function of time. This procedure 
works best when RNA folding is described using two-state approximation. For multidomain 
ribozymes the folding/unfolding kinetics is complex and new tools are required to interpret 
the kinetic data. In a pioneering study Onoa et. al. showed that the rips in FECs 
for the L-21 derivative of Tetrahymena thermophila ribozyme (T. ribozyme), composed of 
multiple domains, are a result of unfolding of individual intact domains that are stabilized 
in the native state by counterion-dependent tertiary interactions. 

The single molecule studies show that the response to mechanical force is a powerful 
tool to analyze the underlying principles of RNA self-assembly. The extraction of unfolding 
pathway using FECs alone is not easy especially when the ribozyme is composed of multiple 
domains 3|. To decipher the unfolding pathways of T. ribozyme Onoa et. al. 2| did a 
series of experiments in which FECs of different independently-folding subdomains were used 
to interpret the order of unfolding of the substructures. In RNA, there is a clear separation 
in the free energies associated with secondary and tertiary interactions. Thus, the FEC 
for a multi-domain ribozyme is, to a first approximation, the union of the FECs for the 
individual domains. Such a strategy can be used to assign a rip of FEC to the unfolding of 
a particular subdomain as long as the contour lengths of two different unfolded motifs are 
not similar. Moreover, it is known that the precise response of RNA to force depends not 
only on the sequence and the native structure but also on how the force is applied [H IE • 
Single molecule experiments can be performed in different modes that includes either force- 
clamp (/ is in constant) 0, 3] or force-ramp (/ varies in a time dependent manner) 
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isl 2o| . Theoretical studies have proposed models for obtaining a number of experimentally 
measurable quantities including FECs for RNA (2ll . |2^ . 2^. Computational studies have 
shown, using RNA hairpin as an example, that the kinetics of unfolding and force-quench 
refolding as well the nature of unfolding depend on the magnitude of / and the loading rate 
(vf) |l7l 2^. These studies show that it is important to complement the single molecule 
studies with computations that can reliably resolve key issues that are difficult to address in 
experiments. 

In this paper, we probe the forced-unfolding dynamics of RNA molecules using a simple 
model. Because these simulations can be used to directly monitor structures in the transition 
from folded to the fully stretched states, unfolding pathways can be unambiguously resolved. 
We introduce the Self- Organized Polymer (SOP) model for RNA that is based only on the 
self-avoiding nature of the RNA and the native structure. We apply the SOP model to 
probe forced-unfolding of a number of RNA structures of varying complexity. Many of the 
subtle features of the variations in the mechanical unfolding as a function of / and Vf can be 
illustrated using P5GA, a simple RNA. For example, we show that the dramatic movements 
in the location of the unfolding transition state occur as rj (or /) is varied. Applications 
to structures of increasing complexity (TAR RNA, prohead RNA from domain Ila of the 
Hapatitus C virus, (j)29 DNA bacteriophage motor) show that discrete intermediates can be 
populated in force-ramp and force-clamp simulations over a certain range of forces. Our 
results show that the response of RNA to force is largely dependent on the architecture of 
the native state. More importantly, we have established that the characterization of the the 
energy landscape requires using force values (or loading rates) over a wide range. 



II. METHODS 

Model: Our goal is to construct a model for obtaining mechanical folding and unfolding 
trajectories for simple RNA hairpins to large ribozymes. The model has to be realistic 
enough to take into account the interactions that stabilize the native fold yet simple enough 
so that the response to a wide range of forces and loading rates can be explored. To this 
end, we introduce a new class of versatile coarse-grained self-organized polymer (SOP) model 
that is particularly well suited for single-molecule force spectroscopy applications of large 
ribozymes and proteins. The SOP model can be used to probe the response of mechanical 
force that is applied by means of force-clamp (constant force), force-ramp, and force-quench. 
The reasons for using SOP model in force spectroscopy applications are the following : 
(1) Forced-unfolding, and force-quench refolding lead to large conformational changes. For 
example, upon application of constant force the end-to-end distance of the RNA changes by 
about (10-100) nm depending on the size of RNA. Currently, single molecule experiments 
(LOT or AFM) cannot resolve structural changes below a few nm. As a result, details 
of the rupture of hydrogen bonds or local tertiary contacts between specific bases cannot 
be discerned from FEC or the dynamics of R alone. Because only large changes in R, the 
variable that is conjugate to force, are monitored it is not crucial to include details of the local 
interactions such as bond-angle and various dihedral angle potentials. (2) We had shown, 
in the context of mechanical unfolding of proteins, that many of the details of unfolding 
pathways can be accurately computed by taking into account interactions that stabilize the 
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native fold [25|. Based on this observation accurate predictions of unfolding forces and 
the location of the unfolding transition states were made for a number of proteins with (3- 



sandwich, a//3, and a-helical folds. Our previous study [25| also suggested that it is crucial 
to take into account chain connectivity and attractive interactions that faithfully reproduce 
the contact map of a fold. (3) Electrostatic interactions are pivotal in RNA. However, under 
physiological condition, counterion concentration is large enough to effectively screen the 
electrostatic repulsion between the phosphate groups. Thus, due to effective screening (small 
Debye length) the repulsive electrostatic potential between phosphate groups is effectively 
short ranged. 

With the above observations in mind we propose, the SOP model for RNA that retains 
chain connectivity and favorable attractive interactions between sites that stabilize the native 
fold. Each interaction center represents the center of mass of a nucleotide. In terms of the 
coordinates {r^, i = 1,2,... A^} of RNA with A^ nucleotide the total potential energy in the 
SOP representation is 

vr = v.ENE + Kir' + vir' 

Af-3 No o 
j=l j=i+3 

N-2 ^ N 

j^^l 2,2+2 .^-^ j=i+3 

The first term is for the chain connectivity. The finite extensible nonlinear elastic (FENE) 
potential [2^ is used with k = 20kcal/{mol-A'^), Rq = 0.2 nm, and rj^j+i is the distance 
between neighboring beads interaction centers i and i + r°^_^_i is the distance in the native 
structure. The use of FENE potential is more advantageous than the standard harmonic 
potential especially when considering forced-stretching because the fluctuations of rj_j+i are 
strictly restricted around r°^_^^ with variation of ±-Ro- The Lennard- Jones potential is used 
to account for interactions that stabilize the native topology. Native contact is defined for 
the pair of interaction centers whose distance is less than Rq = 1.4 nm in the native state 
for \i — j\ > 2. If z and j sites are in contact in the native state. Ay = 1, otherwise A^j = 0. 
We used et = 0.7 kcal/mol for the native pairs, = 1 kcal/mol for non-native pairs. In the 
current version, we have neglected non-native interactions which will not qualitatively affect 
the results because, under tension, such interactions are greatly destabilized. To ensure the 
non-crossing of the chain, we set a = 7A. Only for i,i + 2 pairs we set a* = 3.5 A to 
prevent the flattening of the helical structures when the overall repulsion is large. There 
are five parameters in the SOP force field (fc, Ro, eh, e;, and Rc) [27|. Of these the results 
are sensitive to the precise values of e/j/e/ and Rc- We have discovered that the quantitative 
results are insensitive to Rc as long as it is in the physical range that is determined by the 
RNA contact maps. In principle, the ratio e^/ei can be adjusted to obtain realistic values of 
forces. For simplicity we choose a uniform value of en for all RNA constructs. Surprisingly, 
the SOP force field, with the same set of parameters, can be used to obtain near quantitative 
results for RNA molecules of varying native topology. 
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The time spent to calculate Lennard-Jone forces scales as 0{N'^). Drastic savings in 
computational time can be achieved by truncating the forces due to the Lennard- Jones 
potential for interaction pairs with > {3rfj or 3a), to zero. We refer to the model 
as "Self-Organized Polymer (SOP) model" because it only uses the polymeric nature of 
the biomolecules with the crucial topological constraints that arise by the specific fold. For 
probing forced-unfolding of RNA (or proteins) it is sufficient to include attractive interactions 
only between contacts that stabilize the native state (see Eq. Q). We believe none of the 
results will change qualitatively if this restriction is relaxed i.e., if non-native interactions 
are also taken into account. 

Simulations: Using the SOP model, we simulated the mechanical unfolding and refolding 
of various RNA structures from a simple hairpin to a large ribozyme [N ^ 400). To simulate 
force-ramp experiments we pull a harmonic spring [kg = 28 pN/nm) that is attached to the 
3' end of molecule at a constant speed (f). The time(t)-dependent force acting on the 3' 
end is / = —ks{z — vt) where z is z*'* coordinate of the 3' end. In force-clamp simulations 
a constant force is applied to one end of the molecule while the other end is fixed. Finally, 
in force-quench computations the force on the molecule is reduced to the final value to 
initiate mechanical refolding. In both force-clamp and force-quench setup the dynamics of 
the linker (usually hybrid RNA/DNA handles) is not relevant whereas depending on the 
characteristics of the linkers the dynamics of linker may play an important role in the force- 
ramp experiments [23 |. 

Time scales: Since a typical value for the mass of a nucleotide, m ~ 300g/mol—A00g/mol, 
the average distance between the adjacent nucleotides in the SOP representation of RNA 
is a 5 A, e/i = 0.7 kcal/mol, the natural time is tl = (^2S_)i/2 = 3^5 pg^ Wg ^gg 
tl = 4. Ops to convert simulation times into real times. To estimate the time scale for 
mechanical unfolding dynamics we use a Brownian dynamics algorithm (iil . l29l | for which 
the natural time for the overdamped motion is th = -^^hTi. We used ( = lOOr^^ in the 
overdamped limit, that approximately corresponds to friction constant of a nucleotide in 
water. 

The equations of motion in the overdamped limit are integrated using the Brownian 
dynamics algorithm. The position of a bead i at the time t + h is given by 

x^{t + h)=Xi{t) + ^{Fi{t)+Ti{t)) (2) 

where Fc{t) = — the Newtonian force acting on a bead i, and T(t) is a random force 
on z-th bead that has a white noise spectrum. The autocorrelation function for r{t) in the 
discretized form is 

(r,(t)r,(t + nh)) = ^^5o,n5,i (3) 

where So,n is the Kronecker delta function and n = 0,1,2,.... All the force simulations are 
performed at T = 300 K. For the integration time step h = O.Itl, 10^ integration time 
steps in the over-damped limit {( = lOOr^^) is W^th = W^^^Htl = ^T^s with = 4ps, 
eh ~ 0.7kcal/mol, and ksT ^ O.Qkcal/mol. The system composed of RNA and the spring 
is extended along the force direction by 6x every lO^r^ = 0.47/is integration time steps. We 
chose 6x = 0.003nm so that the pulling speed, v = ^q^j^^ = 6.4 fim/s for h = O.lr^. In 
order to maintain numerical stability neither h nor 6x should be too large. 
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Contact map: In RNA with simple native structures force-induced unfolding pathways can 
be qualitatively predicted from the native structure. To rationalize the simulated unfolding 
pathways it is useful to construct RNA contact maps. We generated the contact maps (Figs. 
I21E1E1) using the definition of native contact. More precisely, we generated a matrix Q where 
matrix elements are 

Q., = e(i?,-r°) (4) 

where r°- is the distance between two nucleotides in the SOP representation of the native 
fold. Representing the native RNA using only the center of mass of each nucleotide as 
the interaction center is a drastic simplification. To ascertain if the SOP representation 
misses any essential feature of the RNA structure we also generated the distance map using 
the heavy atom (C, N, O, P) coordinates. The coarse-grained model captures the important 
interactions on length scales that is larger than about 0.7 nm. For example, the SOP contact 
map (Fig^-C) and the distance map (Fig^) of luud are similar. 

Dynamics of rupture of contacts: The dynamics of RNA unfolding is monitored using 
a number of variables including the time dependence of R and the number of nucleotide- 
dependent native contacts Qi{t) that remain at time t. We define Qi(t) = ^^|j_j|>2) ©(-^C" 
rij{t))Aij where Rc is the cut-off distance for native contacts, rij{t) is the distance between 
i and j-th nucleotide, and Aij = 1 for native contact otherwise Aij = 0. If a certain 
subdomain of the molecule is disrupted and loses its contacts, then the extension of the 
molecule suddenly increases and the mechanical force exerted on the end of the molecule 
drops instantly. These molecular events are reflected as rips in the FEC. When the time 
dependence of the force f{t) or the end-to-end distance R{t) is directly compared with Qi(t) 
using t as a progressive variable to describe unfolding, the direct correlation between sudden 
drops (sudden increase) in the value of f(t) {R(t)) and Qi(t) enables us to unambiguously 
identify the structures involved in the dynamics of rupture of contacts at the nucleotide level. 



III. RESULTS AND DISCUSSION 

Mechanical unfolding of the secondary structural elements of RNA 

The stability of RNA molecule in the native state can be approximated as the sum 
of interactions e*°* = e^'^'^ + Ylk ^k' where ef''^ and e^^*" are interactions that stabilize 
the secondary and tertiary structures, respectively, i refers to the number of secondary 
structural elements, and k labels the tertiary contacts that may be mediated by counterions. 
The contributions from the tertiary interactions are small compared to the energetics 
associated with the secondary interactions (^^ ef''^ ^ ^h'^\ Because of the stability gap 
between the secondary and the tertiary interactions the analysis of FEC for RNA can be 
independently made domain by domain. The hairpin stacks, that can vary in the length 
and sequence, are among the simplest structural motifs. Additional structural complexity in 
RNA arises due to the presence of hairpin loops, bulges, internal loops, internal multi-loops. 
The remarkable structural diversity of RNA secondary structures allows us to probe the se- 
quence and fold-dependent energy landscape using force as a perturbation. Here, we discuss 
the force spectroscopy of relatively simple RNA motifs using four examples. Many aspects 
of the physics of mechanical unfolding of RNA, such as the shifts in the transition state lo- 
cations as r/ is changed, can be understood using these simple structural motifs as examples. 
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Force-induced transitions in a simple hairpin (P5GA): Liphart et. al. showed 
that the P5ab hairpin, the construct in which P5c stem-loop and the A-rich bulge in P5a 
are removed from the P5abc subdomain in T. ribozyme, reversibly folds in an all-or-none 
fashion upon application of constant force. The equilibrium between the native basin of 
attraction (NBA) and the unfolded basin of attraction (UBA) can be shifted by altering 
the value of the constant force, fc- To probe the two-state behavior of hairpins under force 
we used a smaller 22- nt hairpin, P5GA (Protein Data Bank (PDB) id: leor) S^]- For 



the P5GA hairpin simulations over a wide range of forces can be performed in reasonable 
times. The topologically simple hairpin has a single tetra-loop and nine consecutive base 
pairs. In an earlier study we showed, using a minimal three interaction site (TIS) model 
in which each nucleotide is represented by three sites, that the dynamical behavior of P5GA 
under tension is qualitatively similar to P5ab. The much simpler SOP representation of 
P5GA allows us to probe exhaustively the folding and unfolding kinetics of the hairpin that 
is manipulated by force-ramp, force-quench, and force-clarnp. 

Constant force: The hallmark of P5ab [ll| and P5GA [17[, when a constant force is ap- 
plied to either the 3' or the 5' ends, is the observation of bistable kinetics. When a constant 
fc is applied to the 3' end P5GA makes transitions (Fig|21-B) between the UBA [R ^ 8nm) 
to the NBA [R ^ 2?7,m). At fc = 14.0, 15.4, 17.5 pN a large number of transitions occur 
over 45ms duration which suggests that the hairpin dynamics is effectively ergodic. As in 
our previous study the equilibrium constant between the folded and unfolded hairpin 
calculated using a long mechanical unfolding trajectory coincides with an independent en- 
semble average calculation i.e., time averages are roughly equivalent to ensemble averages. 
When fc = 14 pN the residence time in the NBA is much greater than in the UBA while 
at fc = 16.8 pN the UBA is preferentially populated (Fig. 121-B). The population of P5GA 
in the NBA changes when fc is varied can be seen in the histogram (P(i?)) of the end- 
to-end distance R (Fig|2)). At fc = 15.4 pN, which is slightly above the midpoint of the 
NBA UBA transition, several jumps between the NBA and UBA are observed. The 
P{R) distribution reflects the bistable nature of the landscape. The free energy profile with 
respect to R is computed using AF{R) = —ksTlogP^R). From P{R) at fc = 15.4 we can 
obtain the free energy of stability of the folded hairpin with respect to the unfolded state 
using AG ~ fcARuF where ARup is the distance between the folded and unfolded states of 
P5GA. Using ARjjf ~ 6nm we find that AG ~ 13kcal/mol. The Vienna RNA package (3ll ]. 
which uses an entirely different free energy parameters for RNA, gives AG ~ 12.8kcal/'mol. 
This comparison shows that the SOP model can, for simple structures, give accurate results 
for stability. At fc = 15.4 pN, the transition barrier is about ~ l.bkBT. The UBA is 
more populated at this value of fc- The observed transition times are much shorter than 
the residence times in each basin of attraction which is also a reflection of the underlying 
cooperativity of the all-or-none of nature of hopping between UBA and NBA. 

Qualitatively similar results were observed in our previous study using the three inter- 
action site (TIS) model However, the values of the midpoint of the force was about a 
factor of two smaller in the TIS model than in the SOP model. Despite the large differences 
in the nature of the force fields the overall results are robust which suggests that it is the 
underlying native structure that determines the nature of the force-induced transitions in 
simple RNA. Indeed, for simple structures the mechanism of forced unfolding in RNA helices 
are imprinted in the contact map (FigI21-A) which is a two dimensional representation of the 
folded hairpin (Methods). The clustered band in the contact map suggests that P5GA should 
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unfold when a critical number of stacking interactions are unzipped. Thus, the contact map 
for P5GA (and presumably P5ab) is consistent with the observed two-state kinetics. As the 
architecture of the native state becomes more complex it becomes difficult to anticipate the 
unfolding mechanism using the contact map alone (see below). 

Force-ramp: We also performed force-ramp simulations by subjecting the P5GA hairpin 
to a continuously changing force, i.e., varying the loading rate (Methods). The simplicity 
of the SOP model allows us to use values of rj that are comparable to those used in LOT 
experiments. At Vf = 45 pN/s (~ lOr^'^^) the force-extension curves show a transition to the 
UBA at around / ~ 13 pN (Fig|21-A). As the force dynamically increases we observe bistable 
fluctuations in the FEC between the NBA and the UBA just as when force is held constant 
(FigEl-A). The conformational fluctuations between the two states are unambiguously seen 
in the time-dependence of the end-to-end distance {R{t)) (Fig|21-B). As time progresses the 
force is ramped up which results in global unfolding {R ^ 8 nm) for t > 400 ms (FigOl-B). 
During the time scale of simulation we find frequent and sharp transition between the UBA 
and the NBA (Fig. HB). 

The location of the unfolding transition state Ax'^^ for proteins and RNA is often esti- 
mated from force-ramp experiments using the variation of the most probable rupture force 
with Tf ([/*,logrj] plot). The loading rate, t/ = df{t)/dt, and can be accurately estimated 
from the slope of the time dependence of f{t) as a function of time. The slope of f{t) as a 
function of t (FigOl-C) is nearly the same as ^ A;^ x v, where v is the pulling speed. Strictly 
speaking, rj = keff x v with k~^j = + k^^^ + k^^^^^^ and ks, kmoi, and kinker are the 
spring constants of the optical trap, the RNA molecule, and linker, respectively. Typically 
kg <^ kmoi, kiinker, thus keff ~ kg [Sj]- Throughout the paper we obtain the loading rate 



using Tf = kgV. 

From the force distributions, computed at four different loading rates (Fig|21-D), we 
observe that the most probable rupture force (/*) does not increase logarithmically over 
a wide range of loading rates (Fig|21-E). Only if the range of rj is restricted /* changes 
linearly with rj |33|. The location of the transition state (Ax^^) is usually calculated using 
/* «i ( Ax^' ) logrj j33| which may be reasonable as long as rj range is small. However, the 
[/*,logrj] plot is highly nonlinear (Fig. |3]-E). If we use linear regression to analyze the 
[/*,logrj] plot then Ax'^^^ ~ 0.8 nm for the distance between the NBA and the transition 
state. The small value of Ax'§^^ is a consequence of the large variation of Ax]^^ as is 
changed 0, 3^. If the loading rate is varied over a broad range, the r ^--dependence of 



Ax'^^ is manifested as a pronounced convex curvature in the [f\\ogrj] plot (Fig. HE). 
Based on the equilibrium free energy profile F{R) as a function of R [l7[ we expect that 
for P5GA Ax'§^^ ~ 3 nm if ry is small. Indeed, from the constant force simulation results 
in Fig|21-B we find Axp^ ^ Ru-Rf _ xhus, the slope of [/*,logr/] should decrease {Ax'^^ 
increases) as rj decreases. To illustrate the dramatic movement in the transition state we 
have calculate Ax'^^ using /* values for two consecutive values of rj. For example, using 
/* fa 12.5 pN at Vf = 450 pN/s and /* ^ 14.9 pN at rj = 4.5 x 10^ pN/s (a value that 
can be realized in AFM experiments) we obtain Ax'§^^ ~ 4.0 nm. From the values of /* at 
five values of r/ we find that Ax|^'^ can move dramatically (FigOJ-E). In particular, we find 
that Ax^"^ changes by nearly a factor of ten as the loading rate is decreased to values that 
are accessible in LOT experiments (see inset in Fig|21-E). Because of the nearly logarithmic 
variation of Ax'§^^ on Vf over a narrow range of Vf we do not expect Ax^'^ to change 
appreciably if Vf is lowered from AbpN/s to ~ 5pN/s. At high rj or the unfolded state 
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is greatly stabilized compared to the folded state. From Hammond's postulate [35[, gener- 
alized to mechanical unfolding it follows that as fc increases Ax^i"^ should move closer 
to the native state. The simulations are, therefore, in accord with the Hammond's postulate. 

Force-induced transition in hairpins with bulges, internal loops - TAR RNA 
(PDB id : luud) The presence of bulges or internal loops contributes to the bending of 
the stiff helical stack. The enhanced flexibility enables formation of intramolecular tertiary 
contacts or facilitates RNA-protein interactions. Base-sugar, base-phosphate as well as base- 
base contacts are found in these structural elements in the native structure. For example, 
the base group of U23 in HIV-1 TAR RNA (PDB code luud Q, FiggJA), protrudes away 
from the hairpin stack and makes tertiary contacts with z = 38 — 41 (see three-dimensional 
structure in FigB l-B). The inherent flexibility in the bulge region also facilitates interaction 
with ligands ^S^- 

To probe the effect of the U23C24U25 bulge (Fig^-A) on mechanical unfolding we per- 
formed force-ramp and force-clamp simulations. The time dependence of R at fc = 14 
pN, (NBA is the preferred state) shows multiple transitions. Unlike in P5GA some of the 
transitions to UBA involves a small pause in intermediate values of i? (5 — 8) nm which 
is suggestive of a short lived intermediate. Force-clamp simulations also show that R fluc- 
tuates between 2 nm (NBA) and 10 nm (UBA). In addition, there is a signature for an 
intermediate with R between (5 — 8) nm. The free energy profile AF{R) shows that there 
is a high free energy intermediate centered at i? ~ 6 nm that is metastable with respect to 
the NBA and the UBA. At fc = 14 pN the intermediate is less stable than the UBA and 
the NBA. It is likely that the instability of the intermediate state might make it difficult 
for experimental detection. As force is increased to fc (= 14.7 pN), which is very close 
to the critical value at which the stabilities of the folded and stretched states are nearly 
equal, the residence times in the NBA and UBA are similar (see middle panel in Fig. EJ-D). 
The P{R) distribution and AF{R) as a function of R show that fc = 14.7 pN is close to 
the critical value. Interestingly, at fc = 14.7 pN the shallow high free energy intermediate 
is less pronounced than at fc = 14.0 pN (Fig. EJ-D). When force is further increased to 
fc = 15.4 pN the intermediate becomes essentially a part of the UBA. When the hairpins 
are unzipped the presence of bulges and internal loops contributes to the formation of the 
intermediate state when tertiary contacts between these bulges and the rest of the structure 
are disrupted. In the predicted intermediate the first six base pairs (FigEl-A) and additional 
contacts associated with these nucleotides are ruptured (see the representative structures in 
Fig. HD). 

The contact map for TAR RNA has two "clusters" (Fig. EJ-C). The one on the upper right 
corner that is associated with the lower hairpin stack and the other (nucleotides 25 — 31) 
represents the structure from the bulge to the apical end. As the lower stack unfolds the 
force propagation goes along the diagonal for the upper cluster to the lower cluster (FigEl- 
C). We can qualitatively predict the kinetic barriers that oppose forced- unfolding at high Vf 
using the contact map computed from the native structure (Fig. EJ-C). If hairpin unfolding 
occurs in a sequential unzipping manner then we expect (upon applying force to the 3' end) 
the structure to disrupt along the direction specified by the series of red arrows in Fig. 
EJ-C. We predict that upon disruption of the first base pair (G17C45) the other contacts 
involving G17 and C45 ((17,39), (17,40), (17,41), (17,42), (17,43), (17,44), (17,45), (18,45), 
(19,45), and (20,45)) should spontaneously break (upper cluster). Similarly upon rupture 
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of G18-C44 base pair the contacts associated with these nucleotides are disrupted (lower 
cluster). Therefore, the contact map suggests that rupture should occur as force propagates 
along the diagonal direction (red arrows in FigEJ-C). When the contacts associated with the 
lower clusters unravel a high energy intermediate is populated (Fig. IH-C). Only detailed 
simulations can reveal the stability and lifetime of the intermediate. 

The free energy profile as a function of R can also be computed from the contact map using 

AF{R) = -iV,e, X (iV,(T = 300K))/Ng-N, x (kBTlJd) log [ ""ff/Z^y ] {N,{T = 300K)) 

is the average number of contacts at temperature T = 300K at zero force, and Nq is the 
corresponding quality in the PDB structure. The first term is the energetic contribution 
arising from Nc surviving contacts and the second term accounts for the entropy arising from 
the segment in which Nd contacts are disrupted. The energetic contribution using A^^ and 
the entropy is given by the product of Nd and the entropy associated with the freely jointed 
chain like model. We used Kuhn length d ~ 2.0 nm, and the effective nucleotide length of 
single strand RNA Igs ~ 0.59 nm. The chain extension is given hjR = Ro + Ndhs where Rq is 
the end-to-end distance in the folded state. In the presence of constant force the free energy 
profile tilts to UBA. The free energy profile AF{R) at fc ~ 17pN pN is suggestive of an 
intermediate R ^ 6 nm whereas at higher force the signature of the intermediate disappears 
(see also Fig. IH-D). Since we approximated the conformational entropy using freely jointed 
chain model in this exercise, the estimate for the equilibrium critical force does not coincide 
with the SOP simulations which shows that the freely jointed chain model does not estimate 
the entropy of the finite-sized RNA structures. Nevertheless, for simple hairpins the number 
of kinetic barriers or kinetic intermediates can be predicted using this simple analysis based 
on the only knowledge of native contact topology. More recently, Cocco et. al. have 
proposed a similar scheme using the Monte Carlo simulations on the sequence dependent 
free energy profile computed with the Turner's thermodynamic rule js^. A similar analysis 
was previously used to estimate equilibrium unbinding force for proteins |4ol |. 

The FECs obtained using force-ramp simulations at = 4.5 x 10^ pN/s (Fig. El-E) also 
show that one intermediate is present. At this value of rj the presence of an intermediate 
occurs as a rip in the FECs (indicated using an arrow) at / ~ 23 pN. However, when rf 
is lowered by a factor of 10 (see right panel in Fig. I3J-E) there is no signature of a rip at 
/ 15 — 17 pN that corresponds to a pause in a high free energy metastable intermediate. 
As time increases the global unfolding are preceded by fiuctuations between metastable 
intermediate to the folded state, and TAR RNA unfolds in an all-or-none manner at a 
force of about 20 pN. The picture that emerges from force-ramp simulations, namely, the 
presence of an intermediate at high rj and its absence at low rj, is completely consistent 
with constant force simulations. The rj-dependent rupture of RNA structures is a general 
property of self-organized molecules which we explore fully using physical arguments and 
simulations of ribozymes (see below). 

HCV IRES domain II: We consider mechanical unfolding of the 55-nt domain Ila of the 
Hepatitis C Viral (HCV) genome whose NMR structure ji^ is known (PDB code: lp5m). 
The secondary structure map of the hairpin contains bulges and is capped by the UUCG 
tertraloop at the apical end (FigEJA). The domain Ha oligonucleotide adopts a distorted 
L-shaped structure ^| (Fig. El-B) with a relatively fiexible hinge bulge (A53-A57) that is 



stabilized by Mg^+. Just as for the TAR RNA the number of plausible kinetic intermediates 
(in the appropriate force regime) in the NBA UBA transition and the range of force 
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and R values over which they occur can be anticipated using the contact map, which reflects 
the nature of the native fold. The contact map (Fig. El-C) shows that it can be partitioned 
into three distinct clusters that are spatially adjacent. Upon application of force to the 
3' end the rupture of contacts associated with nucleotides Gl and C55 occurs and force 
propagates diagonally (see upper right corner of Fig. 3C). There is a change in the structure 
of contact map with the breaking of contacts involving the base pair (A8-U49) which signals 
the formation of the first intermediate. As a result force propagates along the diagonal 
associated with second cluster. Upon disruption of A8U49 base pair force propagates along 
the diagonal of the sub-contact map (blue line in Fig. |31-C). The second intermediate is 
populated when all the base contacts involving this substructure unravel. Similarly, for lp5m 
we expect population of the third intermediate that opposes forced unfolding associated with 
substructure involving contacts associated with the nucleotides that are in the vicinity of 
the green line in Fig. \^C. 

Explicit force-clamp simulations confirm that there indeed are three intermediates asso- 
ciated with mechanical unfolding of lp5m. We generated ten unfolding trajectories for a 
cumulative 800ms at fc = 21pN. As indicated approximately by the energy profile (Fig. 
inj-C) we find multiple transitions between the structures that are revealed as plateaus in 
R{t) (Fig. El-D). Not all the possible transitions are explicitly observed in all the trajecto- 
ries. This observation may be a reflection of the intrinsic heterogeneity or stochastic nature 
of fluctuations. By averaging the residence time in each basin of attraction over time and 
the initial conditions we obtain P{R) and the associated free energy profile (Fig. |31-D). The 
NBA UBA occurs through a sequence of three intermediates. The barrier separating 
the UBA and the first intermediate {R ~ 5nm) is larger than the subsequent ones. All the 
intermediates, whose populations are fc dependent (see below), are metastable with respect 
to NBA and UBA. The predicted intermediates can be detected experimentally provided 
they are long-lived. If the lifetime of the metastable intermediate is too short then a high 
time resolution would be required. Identical unfolding pathway is also found in force-ramp 
simulations (Fig. El-E) in which we find, in many trajectories, three rips that represent the 
kinetic intermediates. The forces at which these intermediates are populated are higher 
because the loading rate is large. The structures that are populated along the unfolding 
pathway (Fig. EJ-D) are in accord with the predictions based on the secondary structure and 
contact maps. 

Just as found explicitly for TAR RNA the intermediates may not be populated at 
lower fc or rf values. From the energy profile in Fig|Sl-C we find (approximately) that 
the energy associated with a putative intermediate at i? 10 nm (see Fig. El-D) is 
AF{R ^ lOnm) ^ bOksT ^ 200pN ■ nm. Thus, the smallest value of / required to 
transiently populate the R ^ 10 nm intermediate is // ~ AF{R ^ 10nm)/10?T.m ^ 20pN. 
Only if / exceeds 20 pN can these intermediates be significantly populated. At these 
forces the predicted intermediates are metastable with respect to both UBA and NBA. 
Their experimental detection would depend on their life times. Unless relatively high time 
resolution is used in experiments, for practical purposes, mechanical unfolding of domain 
Ila of the HCV IR2S RNA might follow two-state behavior. The simple consideration and 
the simulations might explain the reason for apparent absence of intermediates in the recent 
constant force LOT experiments on TAR RNA performed at / < 15 pN. 

Three-way junction prohead RNA (pRNA) : The next level of complexity in 
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the secondary structural RNA motif is the one that contains an internal multiloop. We 
choose the prohead RNA (PDB code Ifoq ji^) which is part of the 029 DNA packaging 
motor. In the context of the motor pRNA assembles as a pentamer with each monomer 
consisting of two stem- loops that are separated by an internal multiloop (FigEl-A). Even 
at this level of complexity it becomes difficult to predict the kinetic barriers associated 
with force-induced unfolding using the contact map alone. For the 109 nucleotide RNA, 
we generated four mechanical force-ramp unfolding trajectories (FigEl-B). The mechanical 
unfolding of the structure exhibits multiple rips signaling the presence of kinetic barriers 
separating the NBA and UBA (FigjUJ-B). The four different FECs, colored in black, red, 
green, and blue, are distinct. Variations in the rip dynamics from molecule to molecule may 
reflect the heterogeneous nature of the unfolding pathways. Upon force-ramp unfolding of 
the pRNA begins at i? ~ 2 nm in all the trajectories (FigEl-B). The EEC in red has three 
rips at R ~15, 28, and 38 nm while the EEC in green has only one rip at R ~38 nm. For 
pRNA, which is structurally more complex than the simpler three way junction (PSabcAA 
for example 0|) the details of the unraveling mechanism is difficult to extract using the 



FECs alone. To unambiguously extract the unfolding pathways we have calculated the time 
dependence of nucleotide-dependent rupture (Fig. El-C) of individual contacts (Methods). 

Because of the additional stem-loop that branches out of the internal multiloop (EigEl-A), 
there are potentially two routes by which the unfolding of pRNA can proceed. One is 
PI =^ P3 =^ P2, and the other is PI =^ P2 =^ P3. We can rule out the unfolding pathways 
such as P2 =^ PI ^ P3 or P3 =4> P2 =4> PI because of the directional nature of the applied 
mechanical force. It should be stressed that such pathways can emerge if force is applied 
at points other than the ends of the chain. Using the time evolution of the native contact 
of the i*'' nucleotide, Qi{t), we unambiguously pin down the time series of the rupture 
of individual contacts upon force-ramp (FigEl-C). The rupture pattern of the individual 
contacts in all the trajectories is very similar (FigEl-C) despite seeming differences in the 
FECs. From EigEl-C it is clear that the nucleotides in PI unravel {Qi{t)) associated with 
these nucleotides are zero for t > 5ms) early (around t ~ 5ms) which signals the first event 
in pRNA unfolding. Subsequently, Qi{t) for i ~ (70 — 85) are disrupted. This shows that 
the smaller stem-loop in P3 opens, forming a cusp at P ~ 15 nm in the EEC. The bigger 
stem-loop with the large bulge in P2 follows the opening of P3 which is manifested as a rip 
near P ~ 35 nm in the EEC. Thus, the unfolding pathway for pRNA is PI =^ P3 =^ P2. 
Depending on the trajectory, not all the three rips are detected in the FECs although the 
unfolding pathway defined by the subdomains is not altered. The sequence of structural 
changes that accompany the unfolding of pRNA for the trajectory in red in FigEl-B is shown 
in EigEJ-D. Force-ramp simulations of pRNA show that even when the native structure 
is not too complex it is difficult to predict the nature of kinetic barriers without detailed 
dynamical simulations or experiments. It would be interesting to test these predictions 
using LOT experiments. 



IV. CONCLUSIONS 

We have used the self-organized polymer representation of a variety of RNA structures to 
predict their mechanical unfolding trajectories. Constant force and force-ramp simulations 



13 



show that dramatic changes in the force profiles take place as the loading rates and the values 
of the force are varied. If the force is varied over a wide range then regions of the energy 
landscape that cannot be accessed in conventional experiments can be probed. However, in 
order to realize the full utility of the single molecule force spectroscopy it becomes necessary 
to use force in distinct modes (constant force, force-ramp, and other combinations) along 
with reliable computations that can mimic the experimental conditions as closely as possible. 
The simulations of RNA, with diverse native structures, using the SOP model illustrate 
the structural details in the unfolding pathways that are experimentally accessible. It is 
remarkable that the simple native state-based SOP model can quantitatively predict the 
FECs for a number of RNA molecules with varying degree of structural complexity. We 
conclude the paper with the following additional remarks. 

Transition state movements show changes from plastic to brittle behavior: The small size 
and simple architecture of RNA hairpins has allowed us to explore their response to force 
over a range of rj that spans four orders of magnitude. The lowest Vf value is close to those 
used in LOT experiments. A key prediction of our simulations is that the location of the 
transition state for P5GA moves dramatically from about 6 nm at low Vf to about 0.5 nm 
at high Vf (see inset to Fig. 5E). The large value of Aa;-^"^ at low Vf suggests that P5GA 
is plastic while the small Ax'§^^ at high Vf is suggestive of brittle behavior. The mechanical 
properties of RNA structures can be drastically altered by varying the loading rate which is 
reminiscent of the changes in the visco-elastic behavior of polymeric materials that changes 
with frequency. The transformation from plastic to brittle behavior can be captured using 
the fragility index js^] used to describe mechanical unfolding of hairpins. Although we have 
discussed the r/-dependent movement of the transition states using P5GA as an example we 
predict that this result is general and should be observed in other RNA structures as well. 

The predicted movement in the transition state as a function of might help resolve 
the apparent differences in the estimated values of Ax^*^ in proteins using LOT and AFM 
experiments. The typical value of Ax'§^^ for a number of proteins using AFM is on the order 
of 0.5 nm ji^ whereas AxJ"^ for RNase H (the only protein that has been experimentally 
studied using the LOT setup) is about Qnm. We believe that such a large difference is not 
merely due to changes in native topology and stability, which undoubtedly are important. 
Rather, it is due to the variations in Vf. The value of rj in LOT is typically less than 
lOpN/s whereas Vf in AFM varies from (100 - 1000)pA^/s. Our findings here suggest that 
the different loading rates used in the two setups might explain the large differences in the 
values of Ax^'^. Thus, it is important to perform experiments on a given protein using both 
LOT and AFM setup to sort out the loading rate dependence of the location of the transition 
state. 

Forced-unfolding pathways for simple RNA are encoded in the contact map: From the 
contact map it is possible to visualize the directions along which the applied tension prop- 
agates. As illustrated using simple RNA structures the unfolding pathway depends on the 
direction of tension propagation and local structural stability. Using the contact map alone 
one can anticipate the most probable unfolding pathways. Here, we have shown that for 
P5GA, TAR RNA, and HCV IRES domain II it is possible to get a qualitative picture of 
forced-unfolding using the reduced representation of RNA structures in the form of contact 
maps. However, as the structural complexity increases and a number of alternate unfolding 
pathways become possible the simple native structure based method alone is not always 
sufficient in predicting how a particular structure unravels. Such is the case in the unfolding 
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of the three-way junction prohead RNA studied here at constant Vf. 

Limitations of the SOP model: The SOP model is remarkably successful in reproducing 
the unfolding pathways of complex ribozymes in a realistic fashion. Surprisingly, for both 
proteins jl^] and RNA, the SOP model is quite successful in predicting the nature of un- 
folding pathways. The unfolding dynamics of proteins as well as RNA, with size exceeding 
250 residues or nucleotides can be conveniently simulated on a PC in a few days with the 
simulation condition used in this paper {rf ^ 10^ — 10^ pN/s). Quenching the force to zero 
drives the stretched state of molecule close to the native state j27|, which allows us to map 
the folding pathways starting from different initial conditions. The SOP model consisting of 
the polymeric nature and the minimal characteristics of RNA architecture is reasonable in 
visualizing the forced-unfolding and force-quench refolding dynamics of RNA molecules. We 
believe that, when accompanied with the experimental analysis, the SOP model can serve 
as a useful tool that provides insights into the folding/unfolding process of large macro- 
molecules. However, like all models there are certain limitations of the SOP model which 
prevent us from making quantitative predictions of the measurable force-extension curves 
especially for large RNA. The sequence and/or the counterion effect are not explicitly taken 
into account in the SOP model. The neglect of explicit counterions in the simulations, fails 
to capture their specific coordination with RNA, which in turn leads to an underestimate of 
the local stability of the folded structure. Hence in applying the SOP model to investigate 
structures in which counterion-mediated tertiary interactions determine local structures eh 
should be varied. Despite the obvious limitations it is clear that the SOP model is powerful 
enough to provide insights into the structures and pathways that are explored upon appli- 
cation of force. The model can not only be used as a predictive tool (as shown here with 
explicit applications on systems for which experiments are not currently available) but also 
can be used to interpret experimental results. 
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V. FIGURE CAPTIONS 



Figure^: The distance map for TAR RNA (PDB code luud). The distance {vij [A]) 
between all the heavy atoms (C, N, O, P) are computed for the nucleotides (|« — j | > 2 with 
Vij < 2 nm). The scale on the right gives the inverse of distances (l/rij) in color. Small 
rectangular lattices corresponds to the nucleotide unit. The index for the heavy atoms is 
labeled in small font on i and j axis, and the index of the nucleotides {i = 17 — 45) is labeled 
using a large font. Upon coarse graining using the SOP model we obtain the contact map 
in Fig. EJC. 

Figure 121: A. The secondary structure of P5GA hairpin and its contact map. B. The 
time-dependent fluctuations of P5GA hairpin between the folded {R ~ l.bnm) and unfolded 
{R ^ 8nm) states. The end-to-end distance changes spontaneously between two values. The 
force-clamped dynamics of the P5GA hairpin is probed for ~ 45 ms. The histograms P{R) 
at fc = 14.0, 15.4, 16.8 pN are shown in red. The free energy profile AF{R) as a function 
of R for fc = 15.4 pN on the right shows two-state behavior. 

FigureO): Force-ramp unfolding of P5GA hairpin. A. An example of FEC at the loading 
rate = 45 pN/s {ks = 0.07 pN/nm, v = 0.64 fim/s). The data points are recorded every 
50 /is (grey color), but for better illustration the running average is displayed every 500 /is 
(red color). B. The end-to-end distance (R) as a function of time. C. Time dependence of 
/. The loading rate Vf = df /dt is nearly a constant. D. Distribution of unbinding forces 
from 100 trajectories at four loading rates rj = r^, y^^, y'"/' I0r°jr {r°^ = 4.5 x 10'^ pN/s 
[kg = 0.7 pN/nm, v = 6.4 fim/s)). The red star is the rupture force value (/* = 12.5 pN) 
at Vf = 45 pN/s. E. Plot /*, the most probable unfolding force, as a function of Vf. The 
position of transition state AxJ^^ is computed using Ax"^^ = ksT—^jp-. The inset shows 
the variation of Ax'§^^ as a function of rj. 

Figure 0] : Analysis of force-induced transitions in TAR RNA hairpin. A. Secondary 
structure. B. In the three dimensional structure the nucleotides {i = 23, 24) in the bulge 
are in yellow color. C. The contact map and the energy profile at three values of /. The 
plot suggests the presence of one intermediate at i? ~ 6 nm when f ^ 0. D. The results 
of force-clamp simulations at fc = 14 pN, 14.7 pN and 15.4 pN. The dynamics of R is 
probed for 150 ms and the histograms P{R) is given in the shaded color. The corresponding 
free energy profiles computed using AF{R) /ksT = — log P{R) are shown on the right. The 
structures that correspond to the three basins of attraction are displayed at the top. E. 
FECs at the loading rates, rj = 4.5 x 10"^ pN/s {kg = 0.7 pN/nm, v = 6.4 fim/s) (left) and 
r°j: = 4.5 X 10^ pN/s [kg = 0.7 pN/nm, v = 0.64 fim/s) (right) for a number of molecules. 
The arrow in the left panel is the signature of a kinetic intermediate which is absent when 
r/ is reduced. 

Figure El: Force-induced transitions in domain Ila of HCV IRES RNA. A. Secondary 
structure. B. The nucleotides [i = 12 — 16) in bulge in the three dimensional structure 
are in yellow. C. The contact map and the associated energy profile as a function of R 
at three values of fc- The three clustered regions are encircled for clarity. The presence 
of three metastable intermediates at R ~5, 15, 22, 28 nm is indicated by arrows. D. The 
results of force-clamp simulations at fc =18.2, 21.0, and 21.7 pN. The multiple dynamical 
transitions in the hairpin is probed by following the time dependence of R. For each force 
the transitions, which are probed for 800 ms duration in 10 trajectories, are shown for three 
molecules. The distributions P{R) averaged over time and the initial conditions (the dots 
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indicate seven other trajectories) are shown next to the displayed trajectories and the free 
energy profile AF{R) is shown below. The color in the AF{R) profiles are the same as in 
P{R) at each force. The structures of the intermediates, the native conformation, and the 
unfolded structures are explicitly shown. E. FECs at the loading rate, r'j = 4.5 x 10"^ pN/s 
{kg = 0.7 pN/nm, v = 6.4 fim/s) for three trajectories. The arrow indicates the kinetic 
intermediates. 

Figure ini: A. The secondary structure of 109 nt pRNA. B. The four FECs of pRNA at 
rj = 4.5 X 10'^ pN/s {kg = 0.7 pN/nm, v = 6.4 fim/s). C. The rupture history of pRNA 
{Qi(t)) corresponding to FECs in B. Here Qi(t) is the number of contacts that the nucleotide 
i has at time t. The scale on the right represents the number of contacts associated with 
the z*^ nucleotide. Darker shades have larger number of contacts than lighter shades. D. 
Structures of pRNA at each stage of unfolding as time progresses is illustrated. Structures 
in (i) — (iii) are associated with kinetic barriers in the red trajectory in B. 
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